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Abstract 

For a dissipative variant of the two-dimensional Gross-Pitaevskii equation with a parabolic trap 
under rotation, we study a symmetry breaking process that leads to the formation of vortices. The 
first symmetry breaking leads to the formation of many small vortices distributed uniformly near the 
Thomas-Fermi radius. The instability occurs as a result of a linear instability of a vortex-free steady 
state as the rotation is increased above a critical threshold. We focus on the second subsequent symmetry 
breaking, which occurs in the weakly nonlinear regime. At slightly above threshold, we derive a one¬ 
dimensional amplitude equation that describes the slow evolution of the envelope of the initial instability. 

We show that the mechanism responsible for initiating vortex formation is a modulational instability of 
the amplitude equation. We also illustrate the role of dissipation in the symmetry breaking process. All 
analyses are conhrmed by detailed numerical computations. 

Keywords: Nonlinear Schrodinger equation, Bose-Einstein condensates, Vortex nucleation. Dissipative 
Gross-Pitaevskii equation. 
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1 Introduction 

The topic of vortex formation upon rotation of an atomic Bose-Einstein condensate has received a tremendous 
volume of attention during the past 15 years, with many of the relevant results finding their way in main 
archival references on the subject, including the books Eziinni- This is natural, not only because of the 
inherent interest in vortices as fundamental structures in this and more generally in atomic, quantum, 
and superfluids systems |28) . but also because this has been a prototypical way of introducing vortices in 
the system. These studies not only include theoretical works but also numerous experiments, in isotropic 
and anisotropic settings, with few or with many atoms, in oblate or prolate traps in at least four distinct 
experimental groups pioneering the early experiments [iiiiisiin]- Even far more recent experiments, 
relying chiefly on other techniques, including the Kibble-Zurek mechanisms utilize rotation as a way of 
controllably producing vortices of a given (same) charge [SH] . 

ft is then natural to expect a large volume of theoretical literature tackling the relevant theme, ft was 
realized early on that the surface excitations play a crucial role in the relevant “instability” that leads to the 
emergence of vortices [9] . The work of Isoshima and Machida [18] was among the first that recognized the 
complex energetic balance between the different scenarios (e.g. stable, metastable or potentially unstable non¬ 
vortex states, and similarly for vortex bearing states). This metastability opens the potential for hysteretic 
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phenomena, depending on whether, in accordance with the experiment, the rotation frequency was ramped 
up or ramped down, as illustrated, e.g., in Ref. m- Numerous simulations also followed these earlier works, 
including, at different levels, finite temperature considerations. More specifically, both Refs. |26], as well as 
Ref. [31] considered the finite temperature model of the so-called dissipative Gross-Pitaevskii equation (see 
details below), and of ramps therein, as a prototypical system where the nucleation and emergence of vortex 
lattices was spontaneous. On the other hand, the work of Ref. m used the framework of the Hartree-Fock- 
Bogoliubov method (in the so-called Popov approximation) as a means of self-consistently including thermal 
effects, finding that the particular value of the temperature may affect the number of vortices formed. 

From a theoretical perspective, there have been, to the best of our understanding, two distinct schools of 
thought. One of these, based on the work of Ref. [5| (see also importantly the later interpretation of Ref. m, 
for the case of a toroidal trap), is based on computing the Landau criterion threshold, i.e., identifying the 
order of the mode that will be associated with the Landau instability and inferring from that the number of 
vortices that will emerge. A distinct approach pioneered by the work of Stringari and collaborators 
(see also Ref. [38|, as well as the review of the relevant considerations in Ref. |30|) involved the bifurcation 
—from the ground state, be it isotropic or anisotropic— of additional states, beyond the rotation frequency 
that renders neutral (i.e., of vanishing frequency) the quadrupolar mode. These two approaches have both 
been developed in the limit of large chemical potential, yet to the best of our knowledge, they have never 
quite been “reconciled” with each other, aside from a short remark in the work of Ref. |5| suggesting that 
Landau method is more relevant when surface excitations are crucial, while if the instability has a more 
global character (e.g., for smaller atom numbers), then the hydrodynamic approach of Refs. [30113^ 135] 

is more suitable. 

While understanding these two approaches, their similarities and differences, and providing a unified 
perspective of this problem based on them appears to us an intriguing problem for further study, we will not 
pursue it further here. Instead, we will focus on characterizing exactly how a vortex is “born” and migrates 
inwards towards the center of the trap. We will build on our earlier work [8] where we used a multi-scale 
expansion to obtain a reduction of the relevant eigenvalue problem, associated with the vortex forming 
instability. In our case, where the model of choice is the dissipative Gross-Pitaevskii equation (DGPE), we 
argued that there is a true instability, contrary to what is the case with the Hamiltonian case, where the 
eigenvalues simply cross the origin of the spectral plane changing “energy” or “signature” —see the details 
of Ref. [8]. Here, we take this analysis a significant step further, by reducing the relevant dynamics, at the 
periphery of the atomic cloud, to an effective one-dimensional azimuthal strip. 

Remarkably, we find that although the original dynamics pertains to a self-defocusing Gross-Pitaevskii 
equation (GPE), this reduced azimuthal evolution bears a self-focusing character. This trait is manifested 
through the emergence of a modulational instability (MI) against the backdrop of the homogeneous back¬ 
ground. This, in turn, results in a “spike” emerging as subtracted from the background, which finally will 
morph into a vortex initially rotating along the strip and gradually spiraling inwards in accordance with its 
dynamical equation of motion —for vortex motion within the DGPE realm see, e.g., Ref. mQ Our emphasis 
here will be in highlighting the mechanism leading to the vortex formation, offering quantitative comparisons 
of our focusing GPE reduction (and its MI mechanism) with the full two-dimensional numerical results. Our 
presentation is structured as follows. In Sec. [^ we briefly present the mathematical setup. In Sec. [^ we 
discuss the weakly nonlinear analysis and the derivation of the effective one-dimensional self-focusing GPE. 
In Sec.]^ we analyze the MI and compare its predictions to the full system. Finally in Sec.[^ we summarize 
our findings and present some directions for future study. 


2 Mathematical Setup 


Our starting point will be the dissipative variant of the Gross-Pitaevskii equation (DGPE) of the form [33] 
(see also the more recent works of Refs. [331 139] 1 


(7 - i)ut 


1 ^ 

-Au + 



U — \u'\^u — iitrotUe, 


( 1 ) 


^In the latter case, the precession corresponds to an anomalous (negative energy or signature) mode and hence the spiraling 
occurs outwards, but in the presence of rotation this mode becomes normal and similar dynamical equations describe the 
spiraling inwards. 
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where \u{p,0,t)\'^ is the time-dependent two-dimensional density of the atomic condensate cloud within a 
parabolic trap of strength Dtrap and 7 accounts for a phenomenological temperature-dependent dissipation 
effect (see, e.g., Refs. [snEsiisaiiiiiiiiii]). This phenomenological dissipation term provides a prototypical 
way of effectively accounting for the interaction of the condensate with the thermal cloud. Physically relevant 
values of 7 > 0 are of magnitude 1 x 10“^; see, e.g.. Ref. [JT]. Equation 0 is already written in the co¬ 
rotating frame of the trap rotating with frequency firot- The chemical potential ^ is a measure of the strength 
of interaction between atoms, which we assume to be large in comparison to all other parameters in Eq. Q. 
This assumption motivates the following scaling and definitions 

t=-T, p=—^y/%lr, u{p,0,t) = y/flW{r,9,T)-, 

''trap 


where 


— f^rot: 


e = 

2pL 


trap 


< 1, 


so that, in rescaled form, the DGPE may be written as 

(7 -i)WT = e^AW +(l-r^)W- - iDWe; 7 > 0, 


( 2 ) 


where the edge of the atomic cloud, the Thomas-Fermi radius, is now rescaled to r = 1. A radially symmetric, 
vortex-free, steady state W = Wo(r) of Eq. ([^ exists and satisfies 

Worr + -Wor + (1 - r'^)WQ - \Wo\^Wo = 0 , |Wo| ^ 0 as r ^ oo. 
r 

Notice that the steady state prohle is identical to the one of the corresponding Hamiltonian (7 = 0) model. 

Let us use here an approach extending our recent considerations in Ref. [H]. In that work, for increasing 
rotation frequency H, it was observed that the steady state Wq first loses stability to a spatial mode scaling as 
( 7 (g- 2 / 3 )_ instability manifests initially in a large number of small vortices distributed uniformly near 
the Thomas-Fermi radius r = 1. In fact, the relevant surface mode going unstable is the one placing in the 
periphery of the system the number of vortices filling it by spanning their respective healing lengths. This 
behavior was analyzed in Ref. [ 8 ], showing it was due to a linear instability of Wq —within the DPGE set¬ 
ting, although the Hamiltonian solution was still identified as dynamically stable— with increasing rotation 
frequency. Numerical solutions of Eq. ([^ revealed that a subsequent symmetry breaking mechanism causes 
only a fraction of these vortices to persist and be pulled into the bulk of the condensate. For our current 
considerations, for H slightly above threshold, we perform a weakly nonlinear analysis to examine the onset 
of this second symmetry breaking process. While the analysis does not predict the fraction of vortices that 
survive or their eventual fate as they form in the fully nonlinear regime, our analysis accurately captures all 
of the dynamics in the early stages of their development. In particular, we show that the weakly nonlinear 
dynamics of the two-dimensional self-defocusing system ([^ is described by a one-dimensional perturbed 
self-focusing nonlinear Schrodinger equation (NLSE). This is both perhaps intuitively unexpected and at 
the same time crucially relevant to the observed phenomenology. This is because, as our analysis shows, 
the initial pattern selection mechanism responsible for the formation of vortices is a MI of a non-stationary 
uniform solution of the one-dimensional amplitude equation, a mechanism (within the continuum, cubic 
nonlinearity considered herein) restricted to the self-focusing variant of the GPE problem. 


3 Weakly nonlinear analysis and amplitude equations 

Since vortices nucleate near the Thomas-Fermi (r = 1) radius with critical wavenumber m ^ when 

H ~ with toq, H ~ ^iX)^ we rescale Eq. 0 according to 

r = 1 -H 0 = T = W = H 

This way, we are restricting our consideration to the small strip of space near the Thomas-Fermi radius, 
while considering small amplitude solutions (since the density approaches zero near that limit), for longer 
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time scales such that the vorticity is expected to emerge. In these rescaled variables, Eq. (|^ to leading order 
becomes 

(7 -i)wt = Wxx + Wyy - (2x + |wp) w - iVtWy-, (3) 


with 

jw I ^ \/—2x as X —>■ — 00 , and |w| — 0 as x —>■ 00 , 

where w is periodic in y. In arriving at Eq. ([^ from Eq. ([^, the largest terms that have been dropped are 
of order The focus of the analysis and computations herein will be on Eq. 

Writing w = u + iv with u,v G M., we rewrite Eq. (§ as the system 


JUt + Vt 


UxX + Vyy — ( 2 X + + V‘^)U + flVy, 


(4a) 


ivt - Ut 


Vxx + Vyy — ( 2 x + + v"^) V — G,Uy. 


(4b) 


A steady state of Eqs. Q may be written as u = uo(x) and x = 0, where uo(x) is the unique solution of a 
Painleve II equation 

Uq = 2 xuo + Mq, 

with limiting conditions (see, e.g., Ref. m) 


Uq ~ V—2x as X —>■ — 00 , Uq —)■ 0 as X —>■ 00 . 


With respect to the full system (Il)> Uq is the corner layer near r = 1 of the steady state solution lEo(r). 
Next, we let u = Uq{x) + (j){x, y, t) and v = tp{x, y, t) in Eq. 0 to obtain 


+ = 4‘xx + (t>yy- {‘2x + 3ul)(j)-3uo(j)'^ -Uo'ip'^ (j)^+ ^ipy, (5a) 

'yi’t-fpt = ipxx+'^yy- {‘2-x + ul)'ip-2uo4>'il) - -VLcjiy. (5b) 

The steady state of Eq. ([^ is then (/) = '0 = 0. Assuming a perturbation of the form (0,0) = (iA(x), B{x)) 
in Eq. 0 (given the invariance of the solution along the angular variable and the periodicity of the latter, 
we decompose it in Fourier modes) and collecting linear terms, we obtain the eigenvalue problem 


A” — A — {2x + 3 uq)A + rnflBi = Xi^A + B), 


( 6 a) 


B” — B — (2x + Uq)B + mflAi = X{'jB — A). 


( 6 b) 


As in Ref. [8], we set A = 0 in Eq. Q and solve the associated eigenvalue problem for Q{m), yielding 
the neutral stability curve depicted in Fig. 1(a) (see Ref. [ 8 ] for a detailed analysis and full results). We 
denote Uq as the smallest value of 11 at which the steady state loses stability to a perturbation with critical 
wavenumber mo (see Fig. 1(a) I. Then, when 11 = fig + with J ^ 1, numerical computations show that 
3fi(A) ~ 3 (A )/7 ~ O(A^). This is depicted in Fig. 1(b) 


To analyze the slow evolution of this perturbation slightly above threshold, we assume the asymptotic 
expansion 


II = flo + 0 = 501+5^02 + 1^^03, 0 = 501 + 5^-02 + 5^03; 0 < 5 < 1. (7a) 

The expansion in Eq. (|7a[ ) is motivated by the expectation that the bifurcation is of the pitchfork type, for 
which 0, 0 ~ 0{\/n — Uq). The solvability condition is then expected to arise at 0{S^). Recalling that the 
lowest order term omitted from the leading order in Eq. 0 is of order C>(e^/^), we require that 5^ ^ 

We next introduce the slow spatial and temporal scales 

x = XI5, y = Yl6, t = Tl5'^. (7b) 


The spatial scale is motivated by the 0(5) band of wav enum bers that acquires a positive growth rate as 11 
is increased an 0(5^) distance above threshold (see Fig. 1(a) I, while the temporal scale is motivated by the 
corresponding scaling of A in Eq. 0 (see Fig. |l(b)|). Below, we assume that 7 = 0(1) with respect to 5. 
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Figure 1: (a) A depiction of the neutral stability curve D.{m) (solid thick line). As D is increased by 0{S'^) 
above flo, an 0{S) band of wavenumbers around m = toq acquires positive growth rate, (b) The scalings 
of 3fi(A) and 3(A) are depicted as D, is increased by 0(5^) above Qq. The figures here are for illustrative 
purposes only. See Ref. [8] for full results and a detailed analysis. 


Substituting Eq. Q into Eq. ([^, we solve the linear problems at successively higher orders of 6. At 0{6), 
we obtain the linear terms of Eq. (O 


I'Plt+i’lt = (l)lxx + (l}lyy - + 3ul)(l)i Fflo-lpiy, 

'yi’lt-4'lt = i!lxx+i>lyy-{‘^X + ul)4li-D.Qcj)iy. 
We calculate a t-independent solution to Eq. ([^ of the form 


V V'l 


= C{X,Y,T) 


iAi{x) 

Bi{x) 


a^nioy 


+ C.C., 


(8a) 

(8b) 

(9) 


where c.c. denotes the complex conjugate. Here, Ai{x) and Bi{x) are real and satisfy 


1 - 

Bi 


A" — mgAi — {2x + 3uo)Ai + m^Vt^Bi 
B” — rrv^Bi — {2x + u^)Bi + mQ^l^Ai 


= 0, Ai,Bi —>• 0 as a; —)■ ±oo, 


with 


mo « 1.111, Ho « 2.529. 
In addition, we impose the normalization constraint 


/ OO 

A\ + B^dx = 1; Ai,Bi > 0. 

-OO 


(10a) 


(10b) 


( 11 ) 


In Eq. (§, C{X, Y, T) is a complex quantity that describes the slowly modulated envelope of the perturbation, 
while in Eq. (10b), mo is the critical wavenumber that first becomes unstable as H is increased above Hq. 
At we have 

4‘2xx + 4’2yy ~ (2a; + 3Uq) (j)2 + ^0'4’2y = Suo^^j + Uglpi — 2(j)ixX ~ ‘^4‘lyY ~ Ho'i/'iy, (12a) 

1p2xx + i’2yy - {“^x + ul) 1p2 - ^0(t>2y = 2uo(j)i1pi - 2'lpixX - ‘^i’lyY + , (12b) 
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with (j)i and i/'i given in Eq. the terms on the right-hand sides of Eq. ( |l2| ) involve terms proportional 
to C' 2 gi 2 moy^ and [Cp, along with the corresponding complex conjugates. We therefore 

write the solution to Eq. (12) as 


02 

■02 


= 


A22{x) 

-iB22{x) 


J2moy 


Cy 


A2lix) 

-iB2i{x) 


where the equations for ^ 22 ( 3 ;), B 22 (x ),..., are given by 

^22 


B- 


22 


—3mo^i + uoBf 
— ‘ZuqAiBi 


Cx 


ia2iix) 
021 (a;) 


|2 / ^20 


B20{x) 


^ 22 , B 22 0 as X ^ ± 00 , 


(13) 


B 2 


f 2y7io^i — GiqB\ 
I ‘ZttiqBi — 


0^21 

021 


- 2 yi; 

-2B[ 


and 


^20 \ _ f ()UoAl + 2 uoBf 


S' 


20 




0 


, ^215 .^21 —^ 0 X —y iboo, 

0^21,021 —>■ 0 as a; —>■ ±00, 

^20, B20 —>■ 0 as a; —>■ ±00. 


(14a) 

(14b) 

(15) 


In Eqs. (13)-(15l, the linear operator Lmn is defined in Eq. (10a). Since there exists a non-trivial solution 


to the self-adjoint system (10a), for solutions to Eq. (14) to exist, the right-hand sides must each satisfy the 
Fredholm conditions 




(Ai,Bi) ( ] dx = 0. 


(16) 


The second condition in Eq. (16) may be seen from integrating by parts once and applying the boundary 
conditions in Eq. ( |10a ). The first condition may be inferred from the fact that a solution to Eq. (14a) exists 
and is given by {A 2 i,B 2 i) = dmoi^^i, Bi), which may be seen by differentiating Eq. (10a) with respect to 
Too while noting that dG,Q/dmQ = 0. This condition, along with the the normalization constraint in Eq. (11), 
yields the identity 

rAiBidx='^. 

7-00 “0 


With Eq. (16) satisfied, we impose the additional orthogonality constraints 


A^ 


{Ai,Bi)[ ) da; = 0, and 


(7li,i?i) ( 21 ]dx = 0, 


to uniquely specify 7 I 21 , 521, 021 and 02i. The solutions of Eqs. (13)-(15) are depicted in Fig. 
At C7((5^), we have that 


03a:a: + 031/1/ ~ (2a; -I- 3ato)03 -I-flo03y — 

03a;a: + 031/1/ ~ (2a;-I-Wo)03 ~ ^o03y = Ripi 


(17a) 

(17b) 


where and R^ contain the secular terms and S'^(a;) respectively, along with other 

non-resonant terms that we need not consider. For completeness, we give the explicit expressions for the 
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(a) Ai and Bi 


(b) A 22 and B 22 




(c) A 21 and B 21 


(d) 021 and /32i 
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X 

(e) 21.20 and B 20 



Figure 2: Solutions of Eqs. (13)-(15), with Ai, 2 I 22 , ^ 21 , CI 21 and A 20 shown in solid, and Bi, B 22 , B 21 , ^21 
and B 20 in dashed. 
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amplitudes of these secular term: 


dC d^C d^C 

S4, = ("fiAi + Bi) — - i (2021 + ^1) + (2moa2i - ^^0^21 - 2031) + 

d^C 

+ i (O0621 - 2moa2i - Ai) - iil2moBiC+ 

+ i (2^0(3^1020 ~ 3^1022 ~ U0-B1&22) + 3 A^ + AiBi^ |C|^C, 

dC d^C d^C 

= {iBi - iAi) — - ( 2/321 + Bi) -7^ - i ( 2 too / 32 i - 1 ^ 0^21 - + 

d^C 

+ (f2oa2i — 2TO0621 — Bi) ~ ^ 2 'mQAiC+ 

+ [2uo (B1O22 ~ ^1^22 + Bia2o) + ^i-Bi + 3i3i] |C|^C. 

The solution to Eq. 0 may then be written as tp^) = (1^31(0;), B^i{x)) + ..., where A^i and B^i 

satisfy the system 


^31 

Bzi 


zSfp 


(18) 


Applying the Fredholm condition to Eq. (181 


{Ai,B 


we obtain the following amplitude equation for C{X,Y,T): 




dx = 0, 


dC d'^G d'^C d'^C 

(zri - 7 - 2 ) ^ + ^Dxvg^ +DYY^,+aC + = 0 . 


(19) 


The coefficients in Eq. (|19|) are all real, and are given by 

Tl ■ 


2 too „ , 2TOn ^ 

“ > 0, r2 = 1, a = > 0 


a = 


Dxx 

Dxy 


D 


YY — 


rio f2o 

/ oo 

Ai [ 2 mo(— 3 Ai A22 + 3A1A20 — B1B22) + 3 Af + Aii?^] + 

-00 

ill[2uo(—Aii322 + A 20 B 1 + A22ili) + A^B\ + 3i?i] dx > 0, 

/ oo 

Ai[-2a'2i - Ai] + Bi[-2^'i - Bi] dx « 0, 

-00 

/ oo 

Ai[ 2A2 i — 2moQ;2i + fio/32i] + i3i[2i32i ~ 2mo/32i + nott2i] dx « 0, 

-00 

/ oo 

Ai[—Ai — 27710^21 + r2oii2i] + iii[~iii ~ 2moil2i + fioA2i] dx > 0. 

-00 


We note that cr > 0 since we assume that the system is above threshold; i.e., 172 > 0 so that 17 > Ho- With 
the normalization (111 and values of mp and fig given in Eq. (10bI, we numerically obtain the following 
values for the coefficients, accurate to the fifth decimal place: 

Tl « 0.87884, T2 = 1, cr « 0.97671172, a « 0.62184, 

( 20 ) 

Dyy « 0.67615, Dxy 10 ■^^ and Dxx 10 ”®. 

The values above show that, remarkably, Dxx and Dxy are very close to zero. While they may or may 
not be exactly zero, for all practical purposes hereafter, we will indeed set them to 0. In this way, the 
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two-dimensional dynamics in the weakly nonlinear regime of Eq. ([^ reduce to dynamics along only one 
dimension. The same reduction was observed for the one-dimensional dynamics of edge modes in a two- 
dimensional NLSE in the presence of a honeycomb potential [2] • Indeed, this may be an indication (exactly, 
or just approximately so) of an effective “topological protection” [3], a theme of intense recent interest in 
the physics community |35j . The reason we indicate the potentially approximate nature of the topological 
protection for our (toroidal) domain strip is that eventually the ensuing vortices escape inwards towards the 
center of the domain. Nevertheless, exploring this aspect in the context of the present work further is an 
especially appealing aspect for further study. 

Further proceeding with our reduction, by neglecting the Dxx and Dxy terms as indicated above, our 
analysis yields the one-dimensional amplitude equation for the envelope C(Y,T) 


dC d^C 

[iTi - 7) — -I- Dyy F aC + a\cfc = 0. 


( 21 ) 


We make one remark regarding the scaling of 7 with respect to S. Due to the C —>■ —C invariance of Eq. H 


the largest term omitted from the amplitude equation (211 is the quintic term |C'|'^C. This term is of 0(^) 
with respect to the rest of the terms in Eq. (21). Therefore, while we assumed in the analysis that 7 = 0{1) 
with respect to <5, we in fact only require that 7 ^ for Eq. (21) to be valid. Further, as we show in the next 


section, when 7 is small, it is responsible only for the growth of low wavenumber perturbations of a spatially 
uniform state, and the dissipation of high wavenumber perturbations. The symmetry breaking mechanism 
in Eq. (21) that initiates vortex formation in the full system ^ occurs on an 0(1) time scale independent 
of 7. As such, Eq. (21) retains all the orders required to accurately describe the weakly nonlinear dynamics 
of the full system. 

To examine the validity of the weakly nonlinear theory, we solved the two-dimensional system @ numer¬ 
ically on the domain x G [—7.5,22.5], y G [— SOtt/too, SOtt/too] so that exactly 80 wavelengths of the critical 
mode Too fit inside the domain of length Ly = IbOTr/mg. The initial conditions were taken to be of the form 
given in Eq. with an envelope randomly perturbed from unity. That is C(Y,0) = 1 -|- 0.01 * rand(2/), 
where rand(2/) takes on a uniformly distributed random value between 0 and 1 at each discrete point in y. 
The parameters 7 and 6 were taken to be 7 = 0.01 and 6 = 0.04. While realistic values of 7 are typically 
smaller than 0.01, this was purely for demonstration purposes and similar results are obtained for smaller. 


more realistic, values of 7. We also simultaneously solved the one-dimensional amplitude equation (21) on 
the domain Y G [—SOt^Tr/mo, SOiJtt/too]; that is, on a domain of length L = 5Ly, consistent with the scaling 


in Eq. (7b). The comparison of the two sets of results is shown in Fig. Each panel 3(a) -3(f) is arranged 
into a left, center, and right column. In the center column of each figure, we show a surface plot of jwl, 
while in the right column, we show a surface plot of 3(ri;) = tp. Blue (red) regions indicate small (large) 
values in the plotted quantity. In the two plots that make up the leftmost column, we show in red a slice of 
|(/)|/(5 (top) and |'0|/<^ (bottom) taken near a; = 0, corresponding to the vicinity of the Thomas-Fermi radius 
where vortices first form in the original system ([^. Here, (p and ip are the real and imaginary parts of the 
perturbation, respectively, and obey Eq. ([^. In black, we plot the envelope obtained by solving Eq. (21). We 
observe excellent agreement, indicating that the two-dimensional dynamics of Eq. @ in the weakly nonlinear 
regime can indeed be captured by the one-dimensional amplitude Eq. (21). 


The initial symmetry breaking is shown in Fig. 3(b)[ where the uniform (7 = 1 state (see Fig. 3(a)) evolves 
into a spatially periodic state. This is due to a MI in Eq. (21), which will be discussed in Sec. |4| Over a 
relatively shorter time scale, the envelope enters the weakly nonlinear regime, oscillating between a slightly 
localized state (see Fig. 3(c)) and a highly localized state (see Fig. |3(d)| ). This stage can still be accurately 
captured by our effective one-dimensional model. It then enters the fully nonlinear regime (see Fig. 3(e)) 
as two of the localized regions become dominant and results in the formation of two vortices that then get 
pulled into the bulk of the condensate (see Fig. |3(f) ). As seen in Figs. |3(e)] and |3(f)| these vortices manifest 
as dips in the surface of jruj. In this regime, the weakly nonlinear results are no longer applicable. However, 
it is clear that the initial MI (see Fig. |3(b)[ ) in the one-dimensional amplitude Eq. is the symmetry 
breaking mechanism responsible for the initiation of the process leading to the formation of vortices in the 
two-dimensional system (|^. 

To show that Eq. ([2I]) is equivalent to (a dissipative variant of) the self-focusing NLSE when 7 = 0, we 
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(e) t = 8859.375 


(f) t = 9421.875 


Figure 3: Evolution of the perturbation at the periphery of the atomic cloud. Each panel (a)-(f) is arranged 
into a left, center, and right column. The center and right columns depict, respectively, the amplitude and 
imaginary part of the solution, where blue (red) regions indicate small (large) values. The left columns depict 
the comparison of the one-dimensional dynamics of the amplitude Eq. (black) versus the dynamics of 
the full two-dimensional system Eq. ([^ (red). The red is taken from an x-slice of (j) (top subpanel) and ip 
(bottom subpanel), the real and imaginary parts of the perturbation, respectively. Starting from random 
perturbations of the envelope about unity (a), the third mode is selected (b). The pattern oscillates between 
a slightly more localized state (c) and a highly localized state (d), before entering the fully nonlinear regime 
(e) where two vortices are nucleated and pulled into the bulk (f). 
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introduce the rescaled variables 


r = ri(l + 7^)r, Y = ^/D^r|, C{Y,T) = J -B{ri,T); where 7=—, 

\ a Ti 


to obtain 


7 — 1 
1 + 7 ^ 


Br = B^r, + (jB + 2 \BfB. 


Lastly, we multiply Eq. (23) across by 7 + i and scale out the rotation by letting 


5(77, r) = 


( 22 ) 


(23) 


(24) 


to obtain 

At = (b 4” i)Arjr] + d^A + 2(7 + i)|Apj4. (25) 

Setting 7 = 0 in Eq. (251, we see that Eq. is equivalent to the self-focusing nonlinear NLSE. Due to 
rotation {A Ae’'^) and dilation (A —>■ AA, rj —>■ Xrj, and r —^ A^r) invariance the self-focusing NLSE admits 
a family of one-soliton solutions of the form 

As{r], t; v,r) = r sech [r(77 -|- 2 vt) ] 6{r], t) = vrj + (w^ — r‘^)T. 

For 7 ^ 1 in Eq. ( |2^ , a perturbation analysis invoking additional translation {rj —>■ ryg) and Galilean 
{A —>■ and rj ^ rj — 2 ct) symmetries leads to a coupled system of equations for the slow time 

evolution of r{jT) and v{jt) (see, e.g., Refs. [TH [131 [TTl for details). 

Based on the time scales in Figs.[^ we find that, once vortices form, they quickly enter the fully nonlinear 
regime. As such, a detailed analysis of the evolution of a localized soliton solution in the amplitude equation, 
the latter of which is only valid in the weakly nonlinear regime of Eq. (§, is not particularly useful. We are 
presently not aware of a technique (aside from the detailed numerical simulations, such as those of Fig. |3(e)| - 
3(f)[ that could capture this second stage of (large amplitude) symmetry breaking. Instead, we focus on the 
initial symmetry breaking mechanism in Eq. (211 that initiates the formation of vortices in Eq. (|^. As seen 
in Fig. |3(b)] this symmetry breaking does occur in the weakly nonlinear regime of Eq. ([^, and is the result 
of a MI in Eq. (l2l|). We analyze this instability in the following section. 


4 Modulational instability 


In this section we analyze the MI of a spatially homogeneous time-dependent solution of Eq. (21). The 
analysis follows that of Ref. [53] ; see also Ref. [53]. To obtain an exact solution of Eq. (25) without the 
spatial term, we take the ansatz 

A = Ao(T) = /(T)e*®fo), 
where the functions / and g satisfy the ODE’s 

f=^[af + 2f% /(0) = |A(0)|; 

(26) 

g' = 2/^ ff(0) = arg(A(0)). 


The system (26) is solved analytically, yielding 


fir) = —1= _ ; g{T) = -7^iog 

V-2-f Cie-2<^T"^ 27 


— 2 “h C16 


— 2<7-yt 


- err + C2, 


(27) 


where 


Cl = 2 -|- 


|AI(0)| 


2 ’ 


C 2 = 7 ^ log 


27 '=V|Al(0)|2 


-f arg(A(0)). 


A spatially homogeneous solution C'o(T) of Eq. (21) is then given by Eq. (27) and the scalings in Eqs. (22) 

1 


and (|24|. In particular, we calculate that 

2cr 


|Co|^ = 


a 




—2 -I- cqc 




Co = 2 -|- 


2cr 

a|Co(0)|2- 


(28) 
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In what follows, we let C —>■ m, T —>• t, F —>• y, and Dyy —t D for cleaner notation. 

To analyze the stability of uo(t), we introduce the perturbation 

u{y,t) = uo{t){1 +ewit) cos qy)-, e 1. (29) 

Substituting Eq. (29) into Eq. and equating coefhcients of cos qy at the leading order in £, we obtain 

— q^Dw + aw + aluop)!!) + 2'u;] = 0. (30) 


/ • \ *^0 / 
(ZTi — 7) —w + w 


Mo 


Next, noting that uo{t) is a solution of Eq. (211, yields 


{iTi - 7)-^ = -cr - a|uo|^- 

Uq 


(31) 


Substituting Eq. (31) into Eq. ( |30[ ) and simplifying yields 

(jTi — ^)w' — [q"^D — |uo|^] w — a\uo\'^w = 0. 

Now, letting w = Wr + iwi and separating real and imaginary parts, yields the matrix eigenvalue problem 


d 

dt 


1 f 1 [{^ + a)\uo\'^ - q"^D] -Ti [{1 - a)\uo\'^ - q"^D] 

7 ^ + Ti \ti [(1 + a)|MoP - q'^D] 7 [(1 - Q;)|uoP - q'^D] 


The system (32l is non-autonomous due to the time-dependence of |uor given in Eq. (recall Co —'Uo)- 



(32) 


However, since 7 <C 1, we observe by Eq. (26) that |uo(OI evolves on an asymptotically slow time scale as 
long as |moI 7“^/^. The ODE system (32) therefore takes the form w' = M(7 t)w, where M(7<) is the 
two-by-two matrix in Eq. (32) with entries that evolve slowly on an 0(7) time scale. This suggests a WKB 
ansatz for w of the form 

w = v(s) s = jt. (33) 


Substituting Eq. (33) into Eq. (32) and collecting terms at leading order in 7, we find that dr/ds satisfies 


the stationary eigenvalue problem Mv = (dr/ds) v. We may thus identify dr/ds with the eigenvalues of M 
computed with its entries frozen in time. Therefore, w = (w^, Wi)'^ grows (decays) when the eigenvalue of M 
with the largest real part lies on the right (left) half-plane. Scenarios in which the aforementioned eigenvalue 
slowly crosses from the left half-plane into the right half-plane may often lead to the phenomenon of delayed 
bifurcations [ilia. That is, the bifurcation may not become observable until the slowly varying control 
parameter responsible for the eigenvalue crossing is an 0(1) distance past the linear stability threshold. This 
phenomenon is absent here since, as we will observe below, the eigenvalue crossing is not slow. As such, we 


will say that an instability in Eq. (32) has been triggered when the largest eigenvalue acquires zero real part. 


The eigenvalue of the matrix in Eq. (321 with largest real part is given by 


X{q) = 


tr 


tr^ 


— det; 


(34) 


where 


tr = 


27[Kp-g^D] 

72 + tI '■ 


and 


det = + q)KP - q^D] [(1 - a)|Mop - q^D] 


r + 


In the limit of small 7, we see that \{q) is 0 (|moP) and positive when q^ lies in the interval {q/_,q‘/_), with 
= (1 ± a)\uo\‘^/D. This band of positively growing wavenumbers is what is responsible for the symmetry 
breaking mechanism that initiates the formation of vortices. To the left of this band, 3fi(A(g)) is 0(7) and 
positive, while to the right of this band, 3fi(A(g)) is 0{jq‘^) and negative. We thus see that the presence of 
small 7 > 0 is responsible for amplification of the low wavenumbers and dissipation of the high wavenumbers. 
The band of instability, and its O(|uop) positive growth rate, would be present even in the case of 7 = 0. 
However, small 7 still influences pattern formation in Eq. (21) through the growth of |uoP by shifting the 
band of instability towards larger wavenumbers. On a finite domain of length L, in which the shortest 
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Figure 4: (a) The dispersion relation given by Eq. (34| for various values of |uoP (0.25, 1, and 4 from left 
to right) as indicated in the legend. As |uoP increases, the band of unstable wavenumbers broadens, shifts 
to the right, and acquires larger growth rates. The height and shape of the bands depend only weakly on 7 
when 7 is small, (b) The time evolution of |uoP given by Eq. (28). The parameters are 7 = 0.01 and fl 2 = 1. 



O 



(a) amplitude equation vs. ODE 


(b) full PDE system vs. amplitude equation 


Figure 5: (a) Growth rate of the perturbation at the periphery of the atomic cloud. The thick solid blue 
line (using the left axis) depicts |rc(t)| as numerically extracted from the PDE solution of Eq. (21) using the 
perturbation prescribed by Eq. ( [2^ . The dashed red line depicts |w(t)| computed from the linearized ODE 
(32). The light solid line black depicts |uo(t)| (using the right axis). Both the ODE and PDE dynamics show 
that |ui(t)| decays initially until |uo(t)| has increased enough so that the smallest admissible wavenumber 
acquires positive growth rate. This happens approximately when |uo(t)| reaches 0.224 (vertical line). While 
the ODE dynamics seems to exhibit a cumulative phase error with respect to the PDE dynamics, it does 
accurately predict when the instability is triggered, (b) Comparison of ^(w) as extracted from the time 
evolution of the full two-dimensional system (red circles) and that extracted from the same time evolution 
of the amplitude equation (21) (blue solid) as in (a). The horizontal axis in both figures is of slow time. The 


parameters are S = 0.04, 7 = 0.005, SLy « 18.09, and D 2 = 1. 
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admissible wavelength qi = 211 jL may initially lie to the right of the band of instability, positive 7 will 
cause the band to drift rightwards and eventually trigger the instability when = Qi- This occurs when 
|uo| = 0{1). To see why this precludes the delayed bifurcations, we observe that as soon as Ai acquires 
positive growth rate, its growth rate is 0(1} positive. Therefore there is no slow crossing of the eigenvalue, 
and hence no delay. 

It is important to note that the larger the L, where L = SLy, the greater the number of wavelengths of 
the unstable mode(s) the domain can contain. Relating this back to the original system ([^, the farther the 
rotation frequency is set above threshold, the more localized regions form in the weakly nonlinear regime. 
This may lead to more vortices in the fully nonlinear regime being pulled into the bulk. The dependence of 
X{q) on |uoP, along with the time evolution of |uoP, are shown in Figs. 4(a) and 4(b) respectively. 

We illustrate the theory using Fig. where the spatially homogeneous state being perturbed is |uo| = 1. 


The domain length is 5Ly « 18.09 so that, by Fig. 4(a) the third (q « 1.04) and the fourth (q « 1.39) 


modes are the two modes that lie in the band of instability, both having similar growth rates. Consistent 
with the theory, three bumps appear in Fig. |3(b)[ Four bumps may also form given the same parameter set 
and different random initial conditions. 

Lastly, we demonstrate how the growth of |mo| due to positive 7 can intrinsically trigger a MI. On a 
domain of length L « 18.09, we solve Eq. ( [^ with |uo(0)| = 0.005 and 7 = 0.005. The smallest admissible 
wavenumber in this domain is qi = 2 it jL « 0.3473. According to Eq. (34), this wavenumber initially 
lies to the right of the band of instability. Therefore, the spatially homogeneous state is initially stable. 
However, as |uo| increases according to Eq. (28), the band of instability drifts to the right. When |mo| 
increases to approximatel y |uo | « 0.224, 3?(A(gi)) becomes positive, and the mode cos(q\y) begins to grow. 
This is illustrated in Fig. 5(a) where we initialize u as u(y,0} = |uo(0)|(l + 2 x 10“® cos((;i?/)). The figure 


depicts with a thick solid blue line the evolution of |ui(t)| as numerically extracted from the PDE solution of 
Eq. (21) using the perturbation prescribed by Eq. ([2^. The dashed red line depicts the evolution of |w(t)| as 


computed from the linearized ODE (32). The initial (oscillatory) decay is due to <71 lying initially to the right 
of the instability band so that A(gi) lies in the left half-plane with Q(X{qi}} ^ 0. When |uo(t)| (light solid 
black line and right axis) increases to approximately 0.224 (thin vertical line), 3?(A((7i)) becomes positive 
real so that the amplitude of the perturbation begins to increase monotonically, verifying the theory. While 
the ODE prediction appears to exhibit a cumulative error in the phase with respect to the PDE dynamics, 
it does accurately predict when the dynamically instability is triggered. 

We next verify that this intrinsic triggering predicted by the MI analysis is also present in the full two- 
dimensional system ([^. In Fig. |5(b)1 we compare as extracted from the time evolution of the full PDE 
system ([^ (red circles) and that extracted from the same time evolution of the amplitude equation ([^ (blue 


solid) as in Fig. 5(a) The parameters in the full PDE were taken to be <5 = 0.04, 7 = 0.005, 6Ly « 18.09, 


and D2 = 1. We observe excellent agreement between the two dynamics. As in Fig. 5(a), Fig. |5(b)] exhibits 


a slow oscillatory decay followed by fast monotonic growth starting near t = 600. Note that the independent 
variable in the horizontal axis in both figures is the rescaled slow time T = S^t, which has been relabeled t 
in accordance with the notation change T —> t in this section. 

To extract ^(w) from the full PDE data, we first calculate from Eq. © that for a given slice x = Xq, 


’4)i(xo,y) 

2B(xo) 


= Cr(Y, T) cos(moy) - Ci(Y, T) sin(mo2/), 


(35) 


where we have defined Cr = 5i(C'), Ci = 3(C'), and Y = 6y and T = 5‘^t are the slow space and time 
variables, respectively. Because of the separation of scales between y and Y in Eq. (35), the quantity 


ijji(xo,y)/(2B(xo)) to leading order takes the form of a slowly modulated phase-shifted cosine of frequency 
toq) the envelope of which is given by \C{Y,T)\. Next, we calculate from Eq. (29) with u —>■ C,t —?> T, and 
y ^ Y, that 




Here, the evolution of |Co(T)p is given analytically by Eq. ( p^ . Thus, to compute 3fJ(u>), we need only 
calculate numerically the envelope of the quantity tj;i{xo,y)/{2B{xo), which yields \C\. The value of 3fi(i(;) 


may then be extracted by numerically computing the amplitude of the left-hand side of Eq. (36), yielding 
the curve marked by red circles in Fig. |5(b)[ This example shows that the phenomenon, predicted by the 
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MI analysis, of initial decay of a spatial perturbation followed by growth persists not only in the reduced 
amplitude equation, but also in the original two-dimensional PDE system. 


5 Discussion 

In the present work, we have revisited the long studied (not only theoretically and numerically, but impor¬ 
tantly also experimentally) problem of the formation of vortices in the presence of rotation. We have argued 
that while a vast literature exists on the subject, there are still various gaps in our understandings of this 
process, including among other things the weakly (and strongly) nonlinear emergence of a single (or a few) 
vortices that eventually travel inward, settling towards the center of the domain. In order to shed light in the 
weakly nonlinear aspect within this process, we have derived a one-dimensional effective amplitude equation 
as a reduction of a dissipative variant of the self-defocusing two-dimensional GPE with a harmonic trap 
under rotation. Remarkably, this equation turns out to be a self-focusing dissipative variant of the GPE. 
The latter has been shown to undergo modulational instabilities and symmetry breakings that eventually 
result in the formation of solitons that lead to the appearance of the vortices drawn inwards in the original 
(full) problem. This is due to two separate symmetry breaking processes. The first, attributed to a linear 
(modulational) instability of a vortex-free, homogeneous steady state of the dissipative GPE as the rotation 
is increased above a threshold, leads to a large number of “small vortices” nucleating near the edge of the 
condensate cloud. The second, which we can monitor numerically, but which is beyond the realm of our 
weakly nonlinear theory, selects a fraction of these small vortices and pulls them into the bulk of the conden¬ 
sate. Not only were we able to derive an effectively one-dimensional equation describing the weakly nonlinear 
state (its one-dimensionality hinting at an approximate topological insulation of the system’s boundary), but 
we were also able to quantify the modulational instability and illustrate that its temporal and spatial scales 
coincide with the emergence of the pattern formation within the full PDE system. While we could not 
capture the final highly nonlinear step of this destabilization and symmetry breaking process analytically, 
our numerical computations shed considerable light to it. Nevertheless, the latter would be an extremely 
intriguing problem for future study. While the specific pattern selection might be the most difficult step to 
tackle, it would also be interesting to perform an analysis along the lines of Ref. m to derive a system of 
equations of motion for the vortices as they move into the bulk. 

Another key problem worth exploring, as indicated in the introduction, is the reconciliation of the surface 
dynamical picture put forth by Ref. [5] (see also, e.g., for a recent exposition. Ref. [10] and our earlier work 
of Ref. 0) and the bulk hydrodynamic approach of Ref. |34j . Lastly, it would also be relevant to perform an 
analysis of vortex formation in the GPE with an anisotropic potential [23] . In the isotropic case considered 
here, the initial instability leads to a uniform formation of small vortices all around the edge of the condensate 
cloud. In contrast, in the anisotropic case, this uniformity is expected to be broken. An analysis could be 
performed to determine where the first vortices are nucleated, and what the subsequent vortex selection 
mechanism is. These problems are currently under study and relevant progress will be reported in future 
publications. 
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